{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#These are the modules you will need\n",
    "from xlayers import finegrid, layers\n",
    "from xlayers.core import layers_numpy\n",
    "from xlayers.core import layers_xarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import some packages\n",
    "import os\n",
    "os.environ['NUMPY_EXPERIMENTAL_ARRAY_FUNCTION'] = '0'\n",
    "\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import gcsfs\n",
    "\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = 12, 6\n",
    "%config InlineBackend.figure_format = 'retina' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>activity_id</th>\n",
       "      <th>institution_id</th>\n",
       "      <th>source_id</th>\n",
       "      <th>experiment_id</th>\n",
       "      <th>member_id</th>\n",
       "      <th>table_id</th>\n",
       "      <th>variable_id</th>\n",
       "      <th>grid_label</th>\n",
       "      <th>zstore</th>\n",
       "      <th>dcpp_init_year</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AerChemMIP</td>\n",
       "      <td>AS-RCEC</td>\n",
       "      <td>TaiESM1</td>\n",
       "      <td>histSST</td>\n",
       "      <td>r1i1p1f1</td>\n",
       "      <td>AERmon</td>\n",
       "      <td>od550aer</td>\n",
       "      <td>gn</td>\n",
       "      <td>gs://cmip6/AerChemMIP/AS-RCEC/TaiESM1/histSST/...</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AerChemMIP</td>\n",
       "      <td>BCC</td>\n",
       "      <td>BCC-ESM1</td>\n",
       "      <td>histSST</td>\n",
       "      <td>r1i1p1f1</td>\n",
       "      <td>AERmon</td>\n",
       "      <td>mmrbc</td>\n",
       "      <td>gn</td>\n",
       "      <td>gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AerChemMIP</td>\n",
       "      <td>BCC</td>\n",
       "      <td>BCC-ESM1</td>\n",
       "      <td>histSST</td>\n",
       "      <td>r1i1p1f1</td>\n",
       "      <td>AERmon</td>\n",
       "      <td>mmrdust</td>\n",
       "      <td>gn</td>\n",
       "      <td>gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AerChemMIP</td>\n",
       "      <td>BCC</td>\n",
       "      <td>BCC-ESM1</td>\n",
       "      <td>histSST</td>\n",
       "      <td>r1i1p1f1</td>\n",
       "      <td>AERmon</td>\n",
       "      <td>mmroa</td>\n",
       "      <td>gn</td>\n",
       "      <td>gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>AerChemMIP</td>\n",
       "      <td>BCC</td>\n",
       "      <td>BCC-ESM1</td>\n",
       "      <td>histSST</td>\n",
       "      <td>r1i1p1f1</td>\n",
       "      <td>AERmon</td>\n",
       "      <td>mmrso4</td>\n",
       "      <td>gn</td>\n",
       "      <td>gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  activity_id institution_id source_id experiment_id member_id table_id  \\\n",
       "0  AerChemMIP        AS-RCEC   TaiESM1       histSST  r1i1p1f1   AERmon   \n",
       "1  AerChemMIP            BCC  BCC-ESM1       histSST  r1i1p1f1   AERmon   \n",
       "2  AerChemMIP            BCC  BCC-ESM1       histSST  r1i1p1f1   AERmon   \n",
       "3  AerChemMIP            BCC  BCC-ESM1       histSST  r1i1p1f1   AERmon   \n",
       "4  AerChemMIP            BCC  BCC-ESM1       histSST  r1i1p1f1   AERmon   \n",
       "\n",
       "  variable_id grid_label                                             zstore  \\\n",
       "0    od550aer         gn  gs://cmip6/AerChemMIP/AS-RCEC/TaiESM1/histSST/...   \n",
       "1       mmrbc         gn  gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...   \n",
       "2     mmrdust         gn  gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...   \n",
       "3       mmroa         gn  gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...   \n",
       "4      mmrso4         gn  gs://cmip6/AerChemMIP/BCC/BCC-ESM1/histSST/r1i...   \n",
       "\n",
       "   dcpp_init_year  \n",
       "0             NaN  \n",
       "1             NaN  \n",
       "2             NaN  \n",
       "3             NaN  \n",
       "4             NaN  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Look at what CMIP6 data is on the cloud\n",
    "import pandas as pd\n",
    "df = pd.read_csv('https://storage.googleapis.com/cmip6/pangeo-cmip6.csv')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load some of that CMIP6 data\n",
    "df_theta = df[(df.table_id == 'Omon') & (df.variable_id == 'thetao')]\n",
    "uri = df_theta[(df_theta.source_id == 'SAM0-UNICON') &\n",
    "                         (df_theta.experiment_id == 'historical')].zstore.values[0]\n",
    "gcs = gcsfs.GCSFileSystem(token='anon')\n",
    "ds_theta = xr.open_zarr(gcs.get_mapper(uri), consolidated=True)\n",
    "df_salt = df[(df.table_id == 'Omon') & (df.variable_id == 'so')]\n",
    "uri = df_salt[(df_salt.source_id == 'SAM0-UNICON') &\n",
    "                         (df_salt.experiment_id == 'historical')].zstore.values[0]\n",
    "gcs = gcsfs.GCSFileSystem(token='anon')\n",
    "ds_salt = xr.open_zarr(gcs.get_mapper(uri), consolidated=True)\n",
    "df_v = df[(df.table_id == 'Omon') & (df.variable_id == 'vo')]\n",
    "uri = df_v[(df_v.source_id == 'SAM0-UNICON') &\n",
    "                         (df_v.experiment_id == 'historical')].zstore.values[0]\n",
    "gcs = gcsfs.GCSFileSystem(token='anon')\n",
    "ds_v = xr.open_zarr(gcs.get_mapper(uri), consolidated=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find density from temperature and salinity\n",
    "import gsw\n",
    "dens = xr.apply_ufunc(gsw.density.sigma0, ds_salt['so'], ds_theta['thetao'],\n",
    "                      dask='parallelized', output_dtypes=[float, ]\n",
    "                     ).rename('dens').to_dataset()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this function calculates the quantities drF and drC from the depth levels of your model \n",
    "#(here lev_bnds defines the depth of each of the cell boundaries)\n",
    "def finegrid_metrics(levs,lev_bnds):\n",
    "    drF = np.diff(lev_bnds,axis=1)\n",
    "    drC = np.concatenate((np.array([levs[0]]),np.diff(levs,axis=0),np.array([lev_bnds[-1,-1]-levs[-1]])))\n",
    "    return(drF,drC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#actually apply the function\n",
    "fine_drf,fine_drc = finegrid_metrics(ds_theta.lev.values,ds_theta.lev_bnds.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f911ab529e8>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAHwCAYAAAAIIrExAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZxddX34/9dn9mwzWUlIWBKWkAACArIEBcSKVGVRoPJ1w6VUbb+1luq3rWLFqt/qrwtWtFht0a/an9Dit1gVwSVAEFEKiGwhIZCwhoQQsmf2z/ePe2dyZph7596Ze+76ej4e8zh3zvmce9/zycm97/nM57w/IcaIJEmSpPJpqnQAkiRJUqMxCZckSZLKzCRckiRJKjOTcEmSJKnMTMIlSZKkMjMJlyRJksrMJFySJEkqM5NwSZIkqcxMwiVJkqQyMwmXJEmSyswkXJIkSSozk3BJkiSpzFoqHUAaQgjrgU5gQ4VDkSRJUn1bDOyIMS4p5qS6TMKBzilTpsxevnz57EoHIkmSpPq1evVq9u7dW/R59ZqEb1i+fPnse++9t9JxSJIkqY6dcMIJ3HfffRuKPc854ZIkSVKZmYRLkiRJZWYSLkmSJJWZSbgkSZJUZibhkiRJUpmZhEuSJEllZhIuSZIklZlJuCRJklRmJuGSJElSmZmES5IkSWVmEi5JkiSVmUm4JEmSVGYlTcJDCK8JIXwvhLAxhNCT3f4khPDGMdquCCHcFELYGkLYE0J4IITwkRBCcyljkiRJkqpNS6meKIRwBfAZYAvwQ2AjMBd4JXAmcFOi7fnA94Bu4HpgK3AucBVwGnBxqeKSJEmSqk1JkvAQwsVkEvCfAW+NMe4cdbw18bgT+DowAJwZY7wnu/+TwErgohDCJTHG60oRmyRJklRtJj0dJYTQBHwB2AO8fXQCDhBj7Et8exEwD7huKAHPtukGrsh++6HJxiVJkqTaNzAYC/4aHIyVDrdgpRgJXwEsAW4AXgohvAk4msxUk7tjjHeNan9WdnvzGM+1ikwyvyKE0B5j7Mn3wiGEe3McWlZo8JIkSao+O7v7ePe1d/Obp7YVfM57T1vMp849KsWoSqcUSfirsttNwH3AK5IHQwirgItijC9kdx2R3a4d/UQxxv4QwnrgKOAQYHUJ4pMkSVKNufmh54tKwGtNKZLw/bLbDwLrgd8Bfg0cDPw98AbgP8jcnAnQld1uz/F8Q/tnjvfCMcYTxtqfHSE/frzzJUmSVJ227+0b8X1TGP+cQAGNqkQpkvChkoKBzIj3b7PfPxxCeAuZEe8zQginjjE1ZSxDvVc7k3okSZJUUr0Dg8OPP3jGofzF79bXbONS1Al/Kbt9IpGAAxBj3Avckv32pOx2aKS7i7F1jmonSZKkBtPbvy8Jb2uunRHuQpUiCV+T3eaatDOUpE8Z1X7p6IYhhBYyN3n2A0+UIDZJkiTVoL7ESHhbS/0t8l6K6SiryCTNh4cQ2mKMvaOOH53dbshuVwLvAM4Bvjuq7enAVGDVeJVRJEmSVJs2bNnNmk0vq2o9wmObdg0/bm02CX+ZGOOWEML1ZBLrv2JfrW9CCK8nc2PmdvaVJLyBTF3xS0IIVycW6+kAPpttc81k45IkSVL1WbX2Bd597d1FneNIeG6XAycDnwghnA7cTaY6ylvIrIx5WYxxG0CMcUcI4TIyyfhtIYTryCxbfx6Z8oU3kFnKXpIkSXVm5aObiz5n8dxpKURSWSVJwmOMm0MIJ5MZBX8LcAqwE/gR8Dcxxl+Nan9jCOEM4BPAhUAHsI5MMv+lGKOVUSRJkupQT//A8OPl+3dywKwpeVrDSYtnc8bh89IOq+xKNRJOjHErmST68gLb3wm8sVSvL0mSpOrX279vrPV9py3m4hMPrGA0lVN/E2wkSZJUteq96kmhGvcnlyRJUtmNrP/duKloyaajSJIkqXFt29PLt+96kmde2pu33QPP7Ftaph5LDxbKJFySJEmTdvXKdfzrL9YXdU6r01EkSZKkiVvzfP7Fd0ab2tbMcQfOTCma6udIuCRJkiatN3HD5QdOP4QleWp7N4XA6Uvn0TWltRyhVSWTcEmSJE1a8obLs49awAkHz6pgNNXP6SiSJEmatGTpwfYGnutdKHtIkiRJk5YcCW/kqieFcjqKJEmScrp1zWb+5qbVbN3dm7dd8ngjL8JTKJNwSZIk5fQPP1nL2k27ijpnWltzStHUD39NkSRJUk5bdvUU1f6txy9iv86OlKKpH46ES5IkKafkDZe3fOR0Zk9ry9m2rbmJrqmNW3awGCbhkiRJyqknccPlgs4Ok+wScTqKJEmSckqOhHvDZek4Ei5JktSAYoys27yL3b0DeduNLD0Y0g6rYZiES5IkNaDL//23/Odvni24fVOAFut/l4w9KUmS1GD6Bwb5/v2FJ+AA+3dNSSmaxuRIuCRJUoPpHRhkMGYehwDHLOrK2356RwsfOP3QMkTWOEzCJUmSGkxffxx+PKO9he//z1dXMJrG5HQUSZKkBtMzsO9mTCueVIa9LkmS1GD6BvaNhLd5s2VFOB1FkiSpjjy2aSc3P/T8iPreo23f2zf8uNWR8IowCZckSaoTe3sHuPif72Lbnr7xG2c5El4Z9rokSVKdeGrrnqIScIDjDpyZUjTKx5FwSZKkOpFc3XK/Ge284+SD87afO6ON845dmHZYGoNJuCRJUp3oTcwDXzRrCn/yO4dXMBrl43QUSZKkOpEcCW91rndV819HkiSpTiQrorRb9aSqOR1FkiSpysUYueqna/n5o5uJMXe7nT2J0oOOhFc1k3BJkqQqd//T2/jSynVFnTOltTmlaFQK/ookSZJU5TZu7y6qfWtz4PzjrHpSzRwJlyRJqnLJGy7PWDqPj73hiLzt9+/qYM709rTD0iSYhEuSJFW5ZOnBeTPaOXpRVwWjUSk4HUWSJKnKWXqw/jgSLkmSVEEDg5EXdvbkbbN1d+/wY0sP1geTcEmSpAp5fns3F17zS57dtrfgc1qbQ4oRqVz8VUqSJKlCbnpwY1EJOMBcb7isC46ES5IkVcjunv7hx9PampnekT81O3phFxedcEDaYakMTMIlSZIqJFn15A9OP5Q/+Z3DKxiNysnpKJIkSRWSTMJbW5zr3UhMwiVJkiokWXqwzdKDDcXpKJIkSSUWY+TX67eybvOuvO0efm7H8OM2Sw82FJNwSZKkErvx/mf50+t/W9Q5LsLTWPzXliRJKrFfrnux6HOWLZiRQiSqVo6ES5IklVhf4obLFYfOYfHcaXnbrzh0Dq88aFbaYamKmIRLkiSVWLLqydtPPog3H7OwgtGoGjkdRZIkqcR6++PwY+d6ayxeFZIkSSWWHAm36onG4nQUSZKkAj3z0h7+6bbH2bhtb952Dz6zbfix9b81FpNwSZKkAv3DT9fyf+97tqhzHAnXWLwqJEmSCrRhy+6i2s/vbOcVi7pSika1rCQj4SGEDcDBOQ5vijEuGOOcFcAVwClAB7AOuBa4OsY4UIq4JEmSSqlvYN8Nl59885EsmTs1Z9umEDhpyWw6WpvLEZpqTCmno2wHvjjG/pet1xpCOB/4HtANXA9sBc4FrgJOAy4uYVySJEkl0ds/sv738v07KxiNalkpk/BtMcYrx2sUQugEvg4MAGfGGO/J7v8ksBK4KIRwSYzxuhLGJkmSNGnJRXgsPajJqMTVcxEwD7huKAEHiDF2k5meAvChCsQlSZIaWIxx3K+exEh4uzdcahJKORLeHkJ4J3AQsBt4AFg1xvzus7Lbm8d4jlXAHmBFCKE9xtiT7wVDCPfmOLSs8LAlSVKj+497nuYzP3yEHd39BZ/jSLgmo5RJ+ALg26P2rQ8hvDfGeHti3xHZ7drRTxBj7A8hrAeOAg4BVpcwPkmSpDF9+dZ1RSXgzU2Bae3ecKmJK1US/g3gDuBhYCeZBPp/An8A/DiEcGqM8bfZtkN1erbneK6h/TPHe9EY4wlj7c+OkB9fWOiSJKnRbdvTV3DbjtYmLnvNIczoaE0xItW7kiThMcZPj9r1EPDBEMIu4M+AK4G3FPh0YehpSxGbJEnSeJI3XD786Tcwrd31DJWutCczfTW7PT2xb2ikO1fl+s5R7SRJklJl1ROVW9pX2ebsdlpi35rsdunoxiGEFmAJ0A88kW5okiRJMDgYRyzC09oc8rSWSiPtv7Wcmt0mE+qVwDuAc4Dvjmp/OjCVTFWVvJVRJEmSxtM3MMivn9jK3r7ci3H3jxgFD4RgEq70TToJDyEcBWyMMW4dtf9g4MvZb7+TOHQD8AXgkhDC1YnFejqAz2bbXDPZuCRJkt71r7/mV09sHb9hVptTUVQmpRgJvxj4ixDCrcB6MtVRDgXeBHQANwF/N9Q4xrgjhHAZmWT8thDCdWSWrT+PTPnCG8gsZS9JkjRhO7v7ikrAAQ6eM238RlIJlCIJv5VM8vxKMtNPpgHbgF+QqRv+7RjjiEonMcYbQwhnAJ8ALiSTrK8DLge+NLq9JElSsXr7R04zOWPpvLztOztaed+rl6QdlgSUIAnPLsRz+7gNX37encAbJ/v6kiRJY+lNzPWeNbWNf7n0VRWMRhrJiU+SJKku9fXv+8N6W4spj6qLV6QkSapLvQP7KqJ4w6WqjctBSZKkmvPLx7dw04MbGRjMfRtZcil6R8JVbUzCJUlSTXlpdy/v/cZ/05O48XI8roKpauMVKUmSasqGF3cXlYADrDhsTkrRSBPjSLgkSaopySXmD54zlQ+cfmje9gu62jn98PzlCaVyMwmXJEk1JVn/e2HXFN5+8kEVjEaaGKejSJKkmtKXqP/tDZeqVV65kiSppvSMWAnTVEa1yekokiSpKgwMRv7XDQ9wx2MvkLvwIHT37av/3e5IuGqUSbgkSaoKdzz2At+775mizpnW3pxSNFK6/PVRkiRVhS27eotqP2tqK7934oEpRSOly5FwSZJUFZI3XJ537EKueNPyvO27prbS3uJIuGqTSbgkSaoKydKDnVNa2K+zo4LRSOlyOookSaoKI0oPNjvCrfrmSLgkSUrdrp5+nnhhV942T23dM/y4tSWkHZJUUSbhkiQpVY9t2skFX7mT3b0D4zfOarf+t+qcV7gkSUrVjx96vqgEHGD/mVNSikaqDo6ES5KkVO1NLK4zv7Od/Wbkv+Hy2AO7OP+4hWmHJVWUSbgkSUpVX6LqyftfvYQ/OP3QCkYjVQeno0iSpFT1JqqetDrXWwJMwiVJUspGlB5sMfWQwOkokiRpgnr7B/nxQxt54oXdedvd//T24ceOhEsZJuGSJGlCvnXXBj77o9VFndPuSLgEOB1FkiRN0P1PbyuqfQhwzAEzU4pGqi2OhEuSpAlJzvV+0zH7c9i86TnbhgCvOXwuS+ZOK0doUtUzCZckSRPSmyg9+NZXLuJ1y+dXMBqptjgdRZIkTUjfQBx+7A2XUnH8HyNJkiak19KD0oQ5HUWSJI3wyHM7+PufrGHTzu687R7fvK80oSPhUnFMwiVJ0gh//5M1/PzRzUWdM6W1OaVopPrkr62SJGmE57bnHwEf7ehFnSxbMCOlaKT65Ei4JEkaIVl68CtvP56DZk/N2balOXDE/Bk0NYVyhCbVDZNwSZI0QrL04FELO1lsbW+p5JyOIkmSRkiOhLda9URKhSPhkiQ1kB3dfeztHcjbpicxEt5m1RMpFSbhkiQ1iK/cuo6rfrqW/sE4fuMsk3ApHf7PkiSpQfzrL9YXlYBPbWtmSpulB6U0OBIuSVKD2NXTP/x43ox28tUzmdLWzAdOP9SVMKWUmIRLktQAYowjbrj81V++jmbLCkoV46+3kiQ1gIHBSMzORGluCibgUoWZhEuS1AB6k2UHm03ApUpzOookSTVu+54+frZ604jSgqPt7dtXltCKJ1LlmYRLklTDBgcjF371l6zbvKvgc7zZUqo8/xdKklTDtuzqKSoBBzhiwYyUopFUKEfCJUmqYcm53lPbmjn/uIV528+c2sbbTzoo7bAkjcMkXJKkGtabmAe+34x2/uatx1QwGkmFcjqKJEk1rG9g3wqYzvWWaof/WyVJqmHJkfBWq55INcPpKJIkVakbf/MsP3pwI4ODMWeb7Xv7hh87Ei7VDpNwSZKq0LPb9vKn/37/8CqXhbD+t1Q7/N8qSVIVevLF3UUl4ACvP3J+OsFIKrlURsJDCO8CvpX99rIY47+M0WYFcAVwCtABrAOuBa6OMQ6Mbi9JUiNJ3nB55P6dXP76pXnb7z+zg6MWdqUdlqQSKXkSHkI4ELga2AVMz9HmfOB7QDdwPbAVOBe4CjgNuLjUcUmSVEv6Ejdc7t/Vwe84yi3VlZJORwkhBOAbwIvAV3O06QS+DgwAZ8YY3x9j/BhwHHAXcFEI4ZJSxiVJUq1JLsLjDZdS/Sn1/+oPA2cB7wV252hzETAPuC7GeM/QzhhjN5npKQAfKnFckiTVlL4BSw9K9axk01FCCMuBzwP/GGNcFUI4K0fTof03j3FsFbAHWBFCaI8x9pQqPkmSqsHe3gHe983/5q4nXiz4HEfCpfpTkiQ8hNACfBt4Cvj4OM2PyG7Xjj4QY+wPIawHjgIOAVaP87r35ji0bJwYJEmqiJ+t3lRUAg7QNaU1pWgkVUqpRsL/Cngl8OoY495x2g7dur09x/Gh/TNLEZgkSdUkubhOIQ6dN43/cdKBKUUjqVImnYSHEE4iM/r99zHGuyYfEiG7Hbc6aozxhBwx3QscX4JYJEkqqeQy8+9ZsZhPnXvkuOdk6h5IqieTSsIT01DWAp8s8LShke5cxUw7R7WTJKlu9I2qemKCLTWmyd7pMR1YCiwHukMIcegL+FS2zdez+76Y/X5NdvuyVQeySf0SoB94YpKxSZJUdZIj4a3NJuBSo5rsdJQe4F9zHDuezDzxX5BJvIemqqwE3gGcA3x31DmnA1OBVVZGkSTVmme37eXRjTvytlm7edfw47bm5rRDklSlJpWEZ2/C/P2xjoUQriSThP+fUcvW3wB8AbgkhHD1UK3wEEIH8Nlsm2smE5ckSeV29/qtvO1rdxHHvaNpn9YWR8KlRlXyZevHE2PcEUK4jEwyflsI4Toyy9afR6Z84Q1klrKXJKlm/PzRTUUl4ABL5kxLJxhJVa/sSThAjPHGEMIZwCeAC4EOYB1wOfClGIt9G5MkqbKSc70P3286B86emrf98QfN5PVHzk87LElVKrUkPMZ4JXBlnuN3Am9M6/UlSSqnZBL+7lMP5l2nLq5cMJKqnuvgSpJUAqNLD0pSPr5LSJJUAiNLD/rxKim/iswJlySpVuzo7uP///VTPPninrzt7ntq2/BjR8IljcckXJKkPL52+xN8+dZ1RZ3jSLik8fguIUlSHo8+n3/xndHamps44eBZKUUjqV44Ei5JUh69A/uq5r5nxWIOnz89Z9umEHj1YXOZO729HKFJqmEm4ZIk5dHbPzD8+Owj57PisLkVjEZSvXA6iiRJefQlRsK94VJSqfhuIklSHpYelJQGp6NIkhrSLx/fwud+tJotu3rytntxV+/wY0fCJZWKSbgkqSF98aeP8fBzxVU+mdbmx6ak0vBXeklSQ3phnBHw0c45agEHzZmaUjSSGo2/0kuSGlJyrvf//cMVLOyakrNtS3Ow7KCkkjIJlyQ1pN6BfUn4oplTmN/ZUcFoJDUap6NIkhpSXyIJb7PqiaQycyRcklR3NmzZza6e/rxtuvv2LcLTatUTSWVmEi5JqitX3Pgg3/nVU0Wd40i4pHLzXUeSVFf+455nimo/d3obrc0hpWgkaWyOhEuS6sbgYKQnUfXkqIWdedtPa2vhstMPIQSTcEnlZRIuSaobyYonbS1N/OjDr6lgNJKUm9NRJEl1o9eKJ5JqhO9QkqS60dc/ciRckqqV01EkSTXhyRd3c9ODz49Y6XK03b37yhJ6s6WkamYSLkmqer39g1z81bvYvLOn4HMcCZdUzXyHkiRVvY3b9xaVgAMce8DMlKKRpMlzJFySVPWSU1BmT2vjnScflLf97GltXPDKRWmHJUkTZhIuSap6yaon8zs7uPzsIyoYjSRNntNRJElVLzkS3uYNl5LqgEm4JKnq9Q3E4cfecCmpHjgdRZJUUf902zpueXgTMcacbXZ1J0sPmoRLqn0m4ZKkinn0+R38fzevKeqcKa3NKUUjSeXjcIIkqWKe27a3qPbNTYHzjluYUjSSVD6OhEuSKiZ5w+Uph8zmL393ed72+3d1sF9nR9phSVLqTMIlSRXTm7jhcu70do490AV2JDUGp6NIkipmZOlBP5IkNQ5HwiVJqRgcjGzZnX+p+a2J45YelNRITMIlSSW3dXcvF17zS9Zv2V3wOZYelNRIfMeTJJXcLQ8/X1QCDjBneltK0UhS9XEkXJJUcsnFdaa0NjOtPf/HzfL9Z3DJqw5KOyxJqhom4ZKkkusd2HfD5XtOW8yfn7OsgtFIUvVxOookqeSseiJJ+fnOKEkqueRIuFVPJOnlnI4iSSpYjJH7n97Gmud35m330LPbhx87Ei5JL2cSLkkq2C0PP88Hv3NfUee0NoeUopGk2uXwhCSpYHeue7Hoc5YumJFCJJJU2xwJlyQVrC8x1/tVi2dxyNzpedu/aslsTj1kTtphSVLNMQmXJBUsWfXkba86iItOOKCC0UhS7XI6iiSpYFY9kaTS8B1UklSw5HSUNm+4lKQJczqKJIlNO7q55rbHeealvXnb3f/0tuHHjoRL0sSZhEuSuHrlY3znV08VdU5bc3NK0UhS/XMYQ5LE+i27i2o/a2orrzxoZkrRSFL9K8lIeAjhC8CJwFJgLrAXeBK4EfhyjPFlhWVDCCuAK4BTgA5gHXAtcHWMcaAUcUmSCpOsevKxNxzB4fvlLj3YFAInHTKbae3+MVWSJqpU76B/CtwH/BTYDEwjk1xfCfxBCOGUGOPTQ41DCOcD3wO6geuBrcC5wFXAacDFJYpLklSA3oE4/HjFoXN45UGzKhiNJNW/UiXhnTHG7tE7QwifAz4O/CXwh9l9ncDXgQHgzBjjPdn9nwRWAheFEC6JMV5XotgkSeNIjoS3NjtTUZLSVpJ32rES8Kx/z24PT+y7CJgHXDeUgCee44rstx8qRVySpMIkSw+2W/VEklKX9oS+c7PbBxL7zspubx6j/SpgD7AihNAeY+xJMzhJqnc/fOA5rvyvR9i6O//b6eC+2SiOhEtSGZQ0CQ8hfBSYDnSRuVHz1WQS8M8nmh2R3a4dfX6MsT+EsB44CjgEWD3O692b49Cy4iKXpPr05ZXr2LKr8PGMEGBGhzdcSlLaSv1O+1FgfuL7m4H3xBhfSOzrym6353iOof3WvpKkSdq2p6/gtq3NgXefupg509tTjEiSBCVOwmOMCwBCCPOBFWRGwH8TQnhzjPG+Ap9maB3kmLdV5vVOGPMJMiPkxxf4epJUt3oTc73v/vjr8ibYAWhqcil6SSqHVCb+xRg3xRj/EzgbmAN8K3F4aKS762UnZnSOaidJmqC+RNWTjrZmmptCzi8TcEkqn1TvvokxPgk8AhwVQpib3b0mu106un0IoQVYAvQDT6QZmyQ1gp7ESHibN1xKUtUox903C7PboVUwVwLvAM4Bvjuq7enAVGCVlVEkKbeBwcg9G7aypzf/AsPJ0oNWPZGk6jHpJDyEsAzYFmN8ftT+JuAzwH7AL2OML2UP3QB8AbgkhHB1YrGeDuCz2TbXTDYuSapnf/Cte/j5o5sLbj805USSVB1KMRJ+DvC3IYRVwOPAi2QqpJxBpszg88BlQ41jjDtCCJeRScZvCyFcR2bZ+vPIlC+8gcxS9pKkMfT2DxaVgAMcPGdqStFIkiaiFEn4z4CvAacBx5IpLbibTB3wbwNfijFuTZ4QY7wxhHAG8AngQqADWAdcnm0/bmUUSWpUySkmzU2B1xw+N09rmN7ewvtevSTtsCRJRZh0Eh5jfAj4owmcdyfwxsm+viQ1mt5ExZPp7S18870nVTAaSdJEeJeOJNWY5Eh4W4tv45JUi3z3lqQa09Nv2UFJqnXlKFEoSSrQvU++xA9++9yI0e7Rdnb3Dz9ubbbiiSTVIpNwSaoSu3r6ufTau9nV0z9+4yxrf0tSbfLdW5KqxNNb9xSVgAOccsiclKKRJKXJkXBJqhLJKSiLZk7hg2cemrf9vOntnLVsv7TDkiSlwCRckqpEsvTg/M523nXKwRWMRpKUJqejSFKV6LX0oCQ1DN/lJalKJEfCveFSkuqb01EkKWWDg5FPfv8hbn10MzFPu+6+geHH1v+WpPpmEi5JKfvvDVv5t18/VdQ5U9t9e5akeuZQiySl7IVdPUW1n9HewttOPDClaCRJ1cChFklKWbL04NlHzufK847K237W1DamtDWnHZYkqYJMwiUpZckbLrumtLJw5pQKRiNJqgZOR5GklPUO7Lsds9XSg5IkHAmXpEnZ2zvAE1t25W3zzNY9w4+teiJJApNwSZqwp17cw5uvvoMd3f0Fn+MiPJIkcDqKJE3YLQ8/X1QCDrCgsyOlaCRJtcSRcEmaoD29+xbXmTu9nXkz2vO2P2phJxeecEDaYUmSaoBJuCRNULL04KWnHswfv+7wCkYjSaolTkeRpAnqTSThVj2RJBXDTw1JmqBk/e9Wq55IkorgdBRJGqV/YJCfPrKJdZvzlx78zVMvDT+26okkqRgm4ZI0yr/f8wwf/88HizqnrTmkFI0kqR45dCNJo9yXGOEu1CsWzUwhEklSvXIkXJJGSc71fv2R8zli/oy87VccOocjF3amHZYkqY6YhEvSKMnSgxcct4g3HbN/BaORJNUjp6NI0igjq54411uSVHom4ZI0SrL+t1VPJElpcDqKpIaxbvMu/u6WNWzc0Z233eOJ0oRt1v+WJKXAJFxSw7jqZ2u5+eHnizqnvbU5pWgkSY3MIR5JDeO5bXuLan/ovGkce0BXStFIkhqZI+GSGkbyhsu/vegYDttves62zU2BI/fvpMXpKJKkFJiES2oYydKDrzigi2ULrO0tSaoMh3gkNYyRpQd9+5MkVY4j4ZLqwu6efvb2DeRt092XKD1oEi5JqiCTcEk179pfrOfzP350RH3v8Vj/W5JUSX4KSap5/+qoB7IAABwlSURBVPqL9cUl4M1NTG93DEKSVDl+CkmqeTu7+4Yfz5raSlPIvdR8R2sz73v1EqaZhEuSKshPIUk1LzkK/os/P8sEW5JU9ZyOIqnm9Q3E4cfO9ZYk1QI/rSTVtIHByMDgviS8pSn3VBRJkqqFf7OVVLV29fSz8tHNdPfmLj3YN5goO9jSRMgzH1ySpGphEi6pKsUYueRrd/HQszsKPsfa35KkWuEnlqSqtKO7v6gEHGDp/OkpRSNJUmk5Ei6pKiWXmG9raeL8Yxfmbd81pZW3n3xQ2mFJklQSJuGSqlJfouzg7Klt/O3Fx1YwGkmSSsvpKJKqUnIkvLXFmy0lSfXFJFxSVUqOhHvDpSSp3jgdRVLZ/fjBjfzggefoTyyyM9rO7v7hx60m4ZKkOmMSLqmsNu/s5o+/+xv6B3Mn4KO1tzanGJEkSeXn8JKksnrmpb1FJeAAr1u2X0rRSJJUGZMeCQ8hzAHeArwJeAWwCOgFHgS+AXwjxjg4xnkrgCuAU4AOYB1wLXB1jDH38niSalpf4obLw/abzkfPPiJv+/md7Rx34My0w5IkqaxKMR3lYuAaYCNwK/AUMB94K/AvwO+GEC6OMQ4PfYUQzge+B3QD1wNbgXOBq4DTss8pqQ71Jm64nN/ZzjlHL6hgNJIkVUYpkvC1wHnAj5Ij3iGEjwN3AxeSSci/l93fCXwdGADOjDHek93/SWAlcFEI4ZIY43UliE1SlUlWPfGGS0lSo5r0J2CMcWWM8Qejp5zEGJ8Hvpr99szEoYuAecB1Qwl4tn03mekpAB+abFySqlNv/7754JYelCQ1qrSro/Rlt/2JfWdltzeP0X4VsAdYEUJojzH2pBmcpNLp7R/kA9++h1+s20LMc9/lYOJga4tJuCSpMaWWhIcQWoB3Z79NJtxDd2GtHX1OjLE/hLAeOAo4BFg9zmvcm+PQsuKilTRZt699gVvXvFDUOV1TWlOKRpKk6pbmMNTngaOBm2KMtyT2d2W323OcN7TfcghSDdm2p7eo9gfNnsrbTzoopWgkSapuqYyEhxA+DPwZ8CjwrmJPz27HLSQcYzwhx+vfCxxf5OtKmoS+xOqXv3fiAXzuLa/I276lKRBCyNtGkqR6VfIkPITwR8A/Ao8Ar4sxbh3VZGiku4uxdY5qJ6kG9PbvK+/f0dps5RNJkvIo6adkCOEjwJeBh4DXZiukjLYmu106xvktwBIyN3I+UcrYJKUrORJuAi5JUn4lGwkPIfw5mXng9wOvjzFuydF0JfAO4Bzgu6OOnQ5MBVZZGUWqHpt3dPPIxh1526zZtHP4cZtVTyRJyqskSXh2oZ2/Bu4Fzh5jCkrSDcAXgEtCCFcnFuvpAD6bbXNNKeKSNHkPPrOdC/7pTgYGx71NY5gj4ZIk5TfpJDyEcCmZBHwAuAP48Bg3W22IMX4TIMa4I4RwGZlk/LYQwnVklq0/j0z5whvILGUvqQr8bPWmohJwgMVzpqYUjSRJ9aEUI+FLsttm4CM52twOfHPomxjjjSGEM4BPkFnWvgNYB1wOfCnGfEt9SCqnnv59i+EumTuNA2fnT7CPO6CLN75i/7TDkiSppk06CY8xXglcOYHz7gTeONnXl5SuvoF9SfjbTzqIy04/pILRSJJUH5y4KSmv3sRIuDdcSpJUGn6iSsorORLuDZeSJJVGKitmSqp+e3r7ue7up1m/ZXfedndv2FfsyJFwSZJKwyRcalDfuHMDf3vLmvEbJrQ2u8y8JEml4LCW1KDGW3xntOamwImLZ6cUjSRJjcWRcKlB9SVuuPwfJx3I8v07c7YNwIrD5rJo5pQyRCZJUv0zCZcaVG/ihsvXHzmfs5bNr2A0kiQ1FqejSA0qWfWkrbm5gpFIktR4TMKlBpWs/+0Nl5IklZfTUaQ6c++TW/nMD1fzws6evO2Sxy09KElSeZmES3Xmiz97jPuf3lbUOdPafSuQJKmcHP6S6sx4I+CjvebwuRy+3/SUopEkSWNx+EuqM8kbLr/z/pNZPHdqzratzU3M7+woR1iSJCnBJFyqM8nSgwfOnsIBs3In4ZIkqTKcjiLVmb7+OPzYGy4lSapOjoRLNeSZl/aws7s/b5u9fQPDj1ubTcIlSapGJuFSjfjfN63ma6ueKOocR8IlSapOfkJLNeL6/366qPadHS1MaXUlTEmSqpEj4VKN2Nu7b5rJEfNnEPIscjmlrZn3v3qJ01EkSapSJuFSDYgxjqh68uM/eQ1NTS41L0lSrXKYTKoBfQP7Kp60NAUTcEmSapxJuFQDkgvweLOlJEm1z+koUoU9u20vP35wIz39gznbdFt2UJKkumISLlXQwGDkkq/dxdNb9xZ8jiPhkiTVPj/NpQrasqunqAQc4JhFXSlFI0mSysWRcKmCehNTUGZ0tPCuUw7O237W1DbecvyitMOSJEkpMwmXKihZdnDe9Hb+1znLKhiNJEkqF6ejSBWUHAn3hktJkhqHn/pSBVl6UJKkxuR0FCkl37hzPT96YCMDMeZss7unf/hxa7ML8EiS1ChMwqUUbNiym0//4JGizulobU4pGkmSVG38+7eUgue2FVd2MAQ4/7iFKUUjSZKqjSPhUgp6EnO9jztwJp9885F52y/o6mDRzClphyVJkqqESbiUgr5E1ZO509s54eBZFYxGkiRVG6ejSCnoHVH1xBsuJUnSSI6ES0WKMfLSnr68bZLH26z/LUmSRjEJl4qws7uPi665izWbdhZ8jovwSJKk0cwOpCL8fPXmohJwgNnT21KKRpIk1SpHwqUi7OxOTDNpaWJaW/7a3ofPn8E7Tz447bAkSVKNMQmXitA7sG/1y3ecfBCfOveoCkYjSZJqldNRpCL0JkoPesOlJEmaKLMIqQh9I0oP+t9HkiRNjNNRpKyHn9vOI8/tyNvmgWe2DT+26okkSZook3AJuH3tC1x67d1FneNIuCRJmiizCAm4c92Wos9ZOn96CpFIkqRG4Ei4xMgbLo89cCaHzcufYB9/8EzOXLpf2mFJkqQ6ZRIuAb2JGy4vPuEA3nmKtb0lSVJ6nI4iYelBSZJUXmYbEpYelCRJ5eV0FNW1F3f18M+rnuCpF/fkbXf/05YelCRJ5WMSrrr21dsf5+t3rC/qHEfCJUlS2sw2VNeeeGF3Ue2nt7dw0uLZKUUjSZKUUZKR8BDCRcAZwHHAscAM4N9ijO/Mc84K4ArgFKADWAdcC1wdYxwoRVxSsurJh886jOX7d+ZsG0Lg5CWz6ZraWo7QJElSAyvVdJQryCTfu4BngGX5GocQzge+B3QD1wNbgXOBq4DTgItLFJcaXLLqyamHzuXUQ+dUMBpJkqSMUk1H+VNgKdAJfChfwxBCJ/B1YAA4M8b4/hjjx8iMot8FXBRCuKREcanB9Y6oehIqGIkkSdI+JUnCY4y3xhgfizHGAppfBMwDrosx3pN4jm4yI+owTiIvFWpE6cHm5gpGIkmStE8lqqOcld3ePMaxVcAeYEUIoT3G2FO+sFRLfvbIJv7q+w+xZVdv3nbJkfBWR8IlSVKVqEQSfkR2u3b0gRhjfwhhPXAUcAiwOt8ThRDuzXEo75x01b4v37qO57Z3F3VOZ4c3XEqSpOpQiSS8K7vdnuP40P6ZZYhFNWrbnvwj4EnNTYHfO/FAFs6ckmJEkiRJhavGxXqG5gyMO788xnjCmE+QGSE/vpRBqbokq57c+tEzWTizI2fbphBcBVOSJFWVSiThQyPdXTmOd45qJ71M78C+39GmtTfT3uJNl5IkqXZUYnhwTXa7dPSBEEILsAToB54oZ1CqLb39+9ZzanOUW5Ik1ZhKjISvBN4BnAN8d9Sx04GpwCorozSmwcHI/c9sY1d3f9523f3J+t8m4ZIkqbZUIgm/AfgCcEkI4eqhWuEhhA7gs9k211QgLlWBD1/3G374wMaiznG+tyRJqjUlScJDCBcAF2S/XZDdnhpC+Gb28ZYY40cBYow7QgiXkUnGbwshXEdm2frzyJQvvIHMUvZqQDc/9HxR7Rd2ddDSZP1vSZJUW0o1En4ccOmofYdkvwCeBD46dCDGeGMI4QzgE8CFQAewDrgc+FKBK2+qzgwMRvoH9/3Tv+bwuXnbT2tr4dIViwnBJFySJNWWkiThMcYrgSuLPOdO4I2leH3Vh2TZwfaWJr79/pMrGI0kSVJ6nEyrqpFcYt6bLSVJUj0z01HVSI6EW3ZQkiTVs2pcMVN16MFntvNfv312RKI92u7eRO1vR8IlSVIdMwlX6rr7Bnj3tb/mpT19BZ9j2UFJklTPzHSUuo3bu4tKwAFOWjI7pWgkSZIqz5Fwpa4vccPlvBnt/NGZh+ZtP2d6O68/cn7aYUmSJFWMSbhSl5wHPm96O+85bUkFo5EkSao8p6ModZYelCRJGsmMSKmz9KAkSdJITkfRhMUY+d83reaWhzcxGGPOdt19joRLkiQlmYRrwh58djtfv2N9UedMbWtOKRpJkqTa4bCkJmzzjp6i2k9pbeb3TjwwpWgkSZJqhyPhmrBk6cHTl87jcxccnbf9rGltTG/3kpMkSTIj0oQlq550TWnlwNlTKxiNJElS7XA6iibMqieSJEkT40i4xtTTP8CGLXvytnn6pb3Dj9taQtohSZIk1Q2TcL3Mph3dvPEf7+DF3b0Fn+NIuCRJUuHMnPQyP3n4+aIScID5XR0pRSNJklR/HAnXy+zpHRh+PGtqK/NmtOdtv2xBJ2+z9KAkSVLBTML1MsnSg28/+SA+9oZlFYxGkiSp/jgdRS+TrHrS6lxvSZKkkjPD0sv0DsThx20tXiKSJEml5nSUBjIwGLltzWbWbNqZt929T24dfmzVE0mSpNIzCW8gP/jtc3zk+vuLOseRcEmSpNIzw2og9yRGuAt19KKuFCKRJElqbI6EN5DkDZdnLJ3H8v0787Y/eclsjj9oVtphSZIkNRyT8AbSl7jh8rxjF3LhCQdUMBpJkqTG5XSUBpIcCXeutyRJUuWYiTWQ3gHrf0uSJFUDp6PUgade3MPf/WQNz27bm7fdY4nShG0tIe2wJEmSlINJeB24euVj/NdvnyvqnI6W5pSikSRJ0nick1AHnnkp/wj4aItmTuH4g616IkmSVCmOhNeBvsRc789ccDTLF8zI2TaEwCsWdXljpiRJUgWZhNeB5A2Xr1jUxXEHzqxgNJIkSRqPw6F1YETpQaueSJIkVT1Hwqtcd98A3X0Dedv0jKj/bdUTSZKkamcSXsW+e/dTfPoHD9PdNzh+4yzrf0uSJFU/M7Yq9i93PFFUAt7cFOia0ppiRJIkSSoFR8Kr2M7u/uHHMzpaaAq5p5q0tTTxnhWLmTm1rRyhSZIkaRJMwqtYsvTg7R97LbOnmWBLkiTVA6ejVLFk1ZPWZm+4lCRJqhcm4VWsbyAOP3ZxHUmSpPrhdJQK2Ns7wO1rN7O7J3/pweQiPK1NJuGSJEn1wiS8Ai79xt3cvX5rwe1bmgJNTU5HkSRJqhcOr5ZZT/9AUQk4wGH7TU8pGkmSJFWCI+FllrzZsqUpcN6xC/O2n97RwjtOPjjtsCRJklRGJuFllrzZckZHC//wtuMqGI0kSZIqwekoZTay7KDdL0mS1IjMAsssuQCPSbgkSVJjcjpKiTy3bS+f/dEj47bblShL2G7tb0mSpIZkEl4iu3v6uenB54s6x5FwSZKkxmQWWEGvXbZfpUOQJElSBVR0JDyEcADw18A5wBxgI3Aj8OkY40uVjK1YC7o6+Mrbjy+4/bwZ7bxq8awUI5IkSVK1qlgSHkI4FPglsB/wfeBR4CTgT4BzQginxRhfrFR8xZrR0cqbjtm/0mFIkiSpBlRyOso/kUnAPxxjvCDG+BcxxrOAq4AjgM9VMDZJkiQpNRVJwkMIhwBnAxuAr4w6/ClgN/CuEMK0MocmSZIkpa5SI+FnZbc/iTEOJg/EGHcCdwJTgVPKHZgkSZKUtkrNCT8iu12b4/hjZEbKlwI/z/UkIYR7cxxaNvHQJEmSpHRVaiS8K7vdnuP40P6ZZYhFkiRJKqtqXawnZLcxX6MY4wljnpwZIS+8XqAkSZJURpUaCR8a6e7KcbxzVDtJkiSpblQqCV+T3S7Ncfzw7DbXnHFJkiSpZlUqCb81uz07hDAihhDCDOA0YC/wq3IHJkmSJKWtIkl4jPFx4CfAYuCPRh3+NDAN+FaMcXeZQ5MkSZJSV8kbM/+QzLL1XwohvA5YDZwMvJbMNJRPVDA2SZIkKTUVW7Y+Oxp+IvBNMsn3nwGHAl8CTo0xvlip2CRJkqQ0VbREYYzxaeC9lYxBkiRJKreKjYRLkiRJjcokXJIkSSozk3BJkiSpzEKMeVeGr0khhBenTJkye/ny5ZUORZIkSXVs9erV7N27d2uMcU4x59VrEr4e6AQ2lPmll2W3j5b5deuZfVpa9mfp2aelZ5+Wlv1ZevZpadV6fy4GdsQYlxRzUl0m4ZUSQrgXIMZ4QqVjqRf2aWnZn6Vnn5aefVpa9mfp2ael1aj96ZxwSZIkqcxMwiVJkqQyMwmXJEmSyswkXJIkSSozk3BJkiSpzKyOIkmSJJWZI+GSJElSmZmES5IkSWVmEi5JkiSVmUm4JEmSVGYm4ZIkSVKZmYRLkiRJZWYSLkmSJJWZSXgJhBAOCCFcG0J4LoTQE0LYEEL4YghhVqVjq7RsX8QcX8/nOGdFCOGmEMLWEMKeEMIDIYSPhBCa87zOpSGEu0MIu0II20MIt4UQ3pzeT5auEMJFIYSrQwh3hBB2ZPvrO+Ock3q/hRCmhBA+HUJYE0LoDiFsDiH8ewhh+WR+3nIopk9DCIvzXLcxhHBdntdpiD4NIcwJIfx+COE/QwjrQgh7sz/vL0II7w8hjPn54nU6tmL702u0MCGEL4QQfh5CeDrbp1tDCL8JIXwqhDAnxzleozkU059eo+NzsZ5JCiEcCvwS2A/4PvAocBLwWmANcFqM8cXKRVhZIYQNwEzgi2Mc3hVj/LtR7c8Hvgd0A9cDW4FzgSOAG2KMF4/xGn8H/BnwDHAD0AZcAswG/jjG+OVS/TzlEkK4HzgW2EXm51oG/FuM8Z052qfebyGEduDnwGnAPcBK4EDgYqAXOCvG+OtJ/eApKqZPQwiLgfXAb4Ebx3i6h2KMN4xxXsP0aQjhg8A1wEbgVuApYD7wVqCLzPV4cUx8yHid5lZsf3qNFiaE0AvcBzwCbAamAacAJwLPAafEGJ9OtPcazaOY/vQaLUCM0a9JfAG3AJHMhZHc/w/Z/V+tdIwV7p8NwIYC23aS+U/dA5yY2N9B5hedCFwy6pwV2f3rgFmJ/YuBF8m8kS6udD9MoN9eCxwOBODM7M/4nUr2G/CX2XP+A2hK7D8/u//h5P5q+yqyTxdnj3+ziOdvqD4FziKTnDSN2r+ATAIZgQu9TlPrT6/Rwn7mjhz7P5eN/5+8RlPrT6/R8X7eSgdQy1/AIdl/4PVjvHHOIDPithuYVulYK9hHGyg8CX9ftj//zxjHzsoeu33U/m9l9793jHP+Onvs05Xuh0n24ZnkTxhT7zcyieuT2f1LxjhnVfbYayvdXyXq04l8eDR0n46K/ePZ2K9O7PM6LW1/eo1Ork+Pzcb+U6/R1PrTa3ScL+eET85Z2e1PYoyDyQMxxp3AncBUMn+qaWTtIYR3hhA+HkL4kxDCa3PMrRvqz5vHOLYK2AOsyP7pqZBzfjyqTb0qR78dChwErI0xri/wnHqwMITwgey1+4EQwjF52tqn+/Rlt/2JfV6nEzdWfw7xGp2Yc7PbBxL7vEYnbqz+HOI1mkNLpQOocUdkt2tzHH8MOBtYSma+UqNaAHx71L71IYT3xhhvT+zL2Z8xxv4QwnrgKDJ/gVgdQpgGLCIzt3zjGK/7WHa7dFLRV79y9Fsh1/roc+rB67Nfw0IItwGXxhifSuyzT7NCCC3Au7PfJj9IvU4nIE9/DvEaLUAI4aPAdDLz608EXk0mYfx8opnXaIEK7M8hXqM5OBI+OV3Z7fYcx4f2zyxDLNXqG8DryCTi04BXAP9M5s9UPw4hHJtoW2x/2v8Z5ei3RuvrPcBngBOAWdmvM8jcMHcm8PPsB8YQ+3SfzwNHAzfFGG9J7Pc6nZhc/ek1WpyPAp8CPkImYbwZODvG+EKijddo4QrpT6/RcZiEpytkt7GiUVRQjPHTMcaVMcZNMcY9McaHYowfJHPj6hTgyiKebqL92bD9n1WOfquraz3GuDnG+FcxxvtijNuyX6vI/GXr18BhwO9P5KmLaFtzfRpC+DCZqgaPAu8q9vTs1us0K19/eo0WJ8a4IMYYyAwIvZXMaPZvQgjHF/E0XqNZhfSn1+j4TMInZ+g3rK4cxztHtdM+X81uT0/sK7Y/x2s/3m/I9aIc/ea1TuZP0sC/ZL8t5tqt+z4NIfwR8I9kSpe9Nsa4dVQTr9MiFNCfY/IazS87IPSfZBLBOWRuBBziNVqkcfoz1zleo1km4ZOzJrvNNdfo8Ow211ylRrY5u03+KSpnf2bnRS4hc2PSEwAxxt3As8D0EML+Y7xGo/R/OfrNa32foT+3Dl+7jd6nIYSPAF8GHiKTMI61EJfXaYEK7M98vEbHEWN8kswvOEeFEOZmd3uNTlCO/szHaxST8Mm6Nbs9O7x8NbMZZArH7wV+Ve7AasCp2e0TiX0rs9tzxmh/OplKM7+MMfYUeM7vjmpTr8rRb4+TqVW8NISwpMBz6tVQtaMnRu1vyD4NIfw5cBVwP5mEcXOOpl6nBSiiP/PxGi3Mwux2ILv1Gp2c0f2Zj9coWCd8sl+4WE++vjkKmD3G/oPJ3LEcgY8n9neS+e244RfrGfUznsn4i/Wk3m/U8IIIE+jTk4G2Mfafle2bCKxo9D4FPpmN856x/q97naban16j4/fnMmDBGPub2Le4zJ1eo6n1p9foOF8uWz9JYyxbv5rMhfdaMn/+WBEbdNn6EMKVwF+Q+YvBemAnmZqebyLzpnYT8JYYY2/inAvILFPbDVxHZsng88guGQz8Xhx10YYQ/h64nJFL3L6NzPy0Wl22/gLgguy3C4A3kBkxuCO7b0uM8aOj2qfab9nauCvJvEneQ6bs5kHUwtLAFNen2fJZRwG3kekfgGPYV2v2kzHGz47xGg3TpyGES4Fvkhn1upqx51xuiDF+M3GO12kOxfan1+j4stN6/pZMje/HySRx88lU6DgEeB54XYzxkcQ5XqM5FNufXqMFqPRvAfXwBRxIphTfRjL/4E+SuaEm70hGvX+R+Y/5XTJ39m8js+DEC8BPydS9DTnOO41Mgv4Smek8DwJ/CjTnea1Lgf8ms0LpTuB24M2V7oNJ9N2VZH6Dz/W1oRL9RqaizafJ/CWjJ/vv+R/AkZXus1L2KfB+4IdkVnzdlf1ZnwKuB14zzus0RJ8W0J8RuM3rNJ3+9BotqE+PBr5CZmrPFjLzubdnf/YryfEZ7TVamv70Gh3/y5FwSZIkqcy8MVOSJEkqM5NwSZIkqcxMwiVJkqQyMwmXJEmSyswkXJIkSSozk3BJkiSpzEzCJUmSpDIzCZckSZLKzCRckiRJKjOTcEmSJKnMTMIlSZKkMjMJlyRJksrMJFySJEkqM5NwSZIkqcxMwiVJkqQyMwmXJEmSyswkXJIkSSqz/wePBkIlOOPt+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 248,
       "width": 368
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Here we call the first xlayers function, finegrid, that calculates key parameters for rebinning the data\n",
    "drf_finer, mapindex, mapfact, cellindex = finegrid.finegrid(np.squeeze(fine_drf),np.squeeze(fine_drc),[fine_drf.size,10])\n",
    "plt.plot(cellindex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#here we define the new coordinate system\n",
    "thetalayers = np.linspace(-3,31,80)\n",
    "\n",
    "#this is the main function in xlayers: it is being applied to a single column here for demonstration purposes\n",
    "VH = layers.layers_1(ds_v.vo[0,:,100,100].values, ds_theta.thetao[0,:,100,100].values,\n",
    "                     thetalayers, mapfact, mapindex, cellindex, drf_finer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#now we prepare more data to go into the \"layers\" function\n",
    "v_in = ds_v.vo[0,:,:,:]#.transpose('lev','time','j','i')\n",
    "theta_in = dens.dens[0,:,:,:]#.transpose('lev','time','j','i')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "thetalayers = np.linspace(20,30,80)\n",
    "\n",
    "#Here we apply the \"layers\" function to an xarray. 'lev' is the name of the vertical coordinate in our dataset\n",
    "#'Tlev' is the name of the new coordinate\n",
    "v_lay1 = layers_xarray(v_in, theta_in,thetalayers, \n",
    "                           mapfact, mapindex, cellindex, drf_finer, 'lev', 'Tlev')\n",
    "#Note, v_lay1 is thickness weighted (i.e. if the input has units m/s, the output has units m^2/s)\n",
    "#To transform back into non-thickness-weighted coordinates, find the thickness of each layer by using \n",
    "#xr.ones_like(v_in) as input, and divide through by this thickness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: theta_in may not be monotonically ascending/descending\n"
     ]
    }
   ],
   "source": [
    "#Here I basically plot the ROC from the thickness-weighted velocity. \n",
    "#Note: I have not summed with dx because I haven't had time to put it in\n",
    "fig = plt.figure(figsize=(10,7))\n",
    "(v_lay1).cumsum('Tlev').sum('i').plot(x='j',yincrease=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:pangeo1]",
   "language": "python",
   "name": "conda-env-pangeo1-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
